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Abstract 

Citrullus colocynthis is a very drought tolerant species, closely related to watermelon (C. lanatus var. lanatus), an 
economically important cucurbit crop. Drought is a threat to plant growth and development, and the discovery of drought 
inducible genes with various functions is of great importance. We used high throughput mRNA lllumina sequencing 
technology and bioinformatic strategies to analyze the C. colocynthis leaf transcriptome under drought treatment. Leaf 
samples at four different time points (0, 24, 36, or 48 hours of withholding water) were used for RNA extraction and lllumina 
sequencing. qRT-PCR of several drought responsive genes was performed to confirm the accuracy of RNA sequencing. Leaf 
transcriptome analysis provided the first glimpse of the drought responsive transcriptome of this unique cucurbit species. A 
total of 5038 full-length cDNAs were detected, with 2545 genes showing significant changes during drought stress. 
Principle component analysis indicated that drought was the major contributing factor regulating transcriptome changes. 
Up regulation of many transcription factors, stress signaling factors, detoxification genes, and genes involved in 
phytohormone signaling and citrulline metabolism occurred under the water deficit conditions. The C. colocynthis 
transcriptome data highlight the activation of a large set of drought related genes in this species, thus providing a valuable 
resource for future functional analysis of candidate genes in defense of drought stress. 
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Introduction 

Water is essential for plant growth in modern agriculture [1]. 
Drought delays the development of crops, and strongly affects 
morphology, as well as physiological processes like transpiration, 
photosynthesis, respiration and translocation of assimilates [2]. 
Drought avoidance can be achieved through morphological 
changes in plants, such as decreased stomatal conductance, 
reduced leaf area, and extensive root systems [3]. Drought 
tolerance is achieved by physiological and molecular mechanisms, 
including osmotic adjustment, antioxidant and scavenger com- 
pounds [4] . Both strategies involve the induction of specific genes 
and proteins, such as dehydrins (dehydration-induced proteins), 
key enzymes for osmolyte biosynthesis, and detoxification enzymes 
[5,6]. 

Plant species have developed diverse strategies to adapt and 
thrive in all kinds of climates and terrains to deal with extreme 
changes in the environment. These strategies are supported by 
rich and complex metabolic networks that enable the plant to 
synthesize a wide range of compounds. Plant responses to abiotic 
stresses involve interactions and crosstalk between many molecular 
pathways. High throughput screening techniques such as tran- 
scriptome sequencing have been used to study the adaptability of 
plants to drought [7] . This led to the discovery of many drought 
related genes. For example, PIP aquaporins were found to fine- 



tune the environment in response to declining water availability 
[8]. However, few natural allelic variants have been cloned for 
drought related traits, so QTL, RNA sequencing, and other trait 
isolation methods are needed to improve methodology for 
exploring complex multivariate data [9, 1 0] . 

The cucurbit family is a large family with several economically 
important species, such as watermelon {Citrullus lanatus), melon 
{Cucumis melo), cucumber [Cucumis sativus) and several Cucurbita 
species with edible fruits [1 1] . Citrullus colocynthis (L.) Schrad 
(2n = 2x = 22), the bitter apple, closely related to domesticated 
watermelon [Citrullus lanatus var. lanatus), is a very drought- 
tolerant perennial herbaceous species in the Cucurbitaceae family 
[12]. It can survive arid environments by maintaining its water 
content under severe stress conditions. C. colocynthis is an 
important medicinal plant and a source of valuable oil [1 2] . Its 
seeds were found in several early Egyptian, Libyan and Near 
Eastern sites from about 4000 BC. This species grows in sandy 
areas throughout northern Africa, southwestern Asia and the 
Mediterranean region [11,12]. The species has been used as a 
model to elucidate the function of genes implicated in the stress 
response ultimately leading to enhancement of stress tolerance in 
cucurbit crops through genetic manipulation. Si et al. [13] found 
dynamic gene expression changes in C. colocynthis root tissues 
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using cDNA amplified fragment lengtli polymorpliism (cDNA- 
AFLP) teclinique. 

Several research groups liave used next generation sequencing 
teclinologies to study gene expression profiles in species of the 
cucurbit family. Guo et al [14] used 454 sequencing technology to 
study the comprehensive profUe for watermelon fruit flesh tissues, 
Grassi et al [15] studied carotenoid pathway regulators in ripening 
watermelon fruit. The draft genome of watermelon (C. lanatus, 
2n = 2x = 22, —425 Mb) was analyzed by Guo et al [16] using 
three different watermelon subspecies. Comparative genomic 
analysis provided an evolutionary scenario for the origin of the 
1 1 watermelon chromosomes derived from a 7-chromosome 
paleohexaploid eudicot ancestor. The genome sequence of 
cucumber {Cucumis sativus, 2n = 2x=14) has been completed, 
and the genome of melon (Cucumis melo, 2n = 2x = 24) is being 
sequenced under the Spanish Genomics Initiative (MELONO- 
MICS) [17,18]. Liu et al [19] used sequencing techniques to 
identify conserved and novel miRNA in watermelon, while 
Wincker [20] used comparative analysis of genomes between 
watermelon and sweet orange to detect the traits related to their 
domestication. 

Here, high-throughput sequencing of the leaf transcriptome 
from C. colocynthis provides a glimpse at drought related genes in 
this uniquely drought tolerant cucurbit species. The lack of 
extensive genomic and functional genomic resources in Citrullus 
species has hampered research and breeding of the cultivated and 
economically important C. lanatus. This study should facilitate the 
identification of valuable multiple genes, needed for complex 
interactions of plant species with the environment. 

Materials and Methods 

Plant Materials and RNA Extraction 

C. colocynthis seedlings were grown in Sunshine Mix 7^8 under 
a 16 h light/8 h dark photoperiod at 26°C day, 22°C night 
temperature. Seedlings with 2-3 true leaves (2-3 week old) in 
50 ml containers were exposed to drought by withholding water 
for 0 (day 1, D 1), 24 (day 2, D2), 36 (day 3, D3) or 48 hours (day 4, 
D4). True leaf samples were collected each day at noon, flash 
frozen in liquid nitrogen and stored at — 80°C. RNA was 
subsequentiy extracted using the TRIzol method [21]. 

Preparation of cDNA Library and Sequencing 

lUumina sequencing was performed at the HudsonAlpha 
Institute of Biotechnology (Huntsville, AL) following manufactur- 
er's instructions. RNA-Seq reads were first processed to remove 
rRNA sequence contamination. First strand cDNA was synthe- 
sized with reverse transcriptase and random primers using the 
small fragment RNAs as template. Second strand cDNA was then 
synthesized followed by phosphorylation by T4 DNA polymerase. 
The cDNA fragments were 3' adenylated and ligated to the 
lUumina's paired-end adapters. Enrichment of cDNA templates 
was conducted following fifteen cycles of PGR amplification. In 
total, over 20*4Mp were sequenced for mapping assembly and 
differential expression analysis. Raw sequence data are available 
for download at NCBI Sequence Read Archive (SAMN02769576, 
SAMN02769577, SAMN02769578, and SAMN02769579). 

Assembly 

Raw sequencing data were filtered Z (0.05) using the CLC 
Genomics workbench (CLC Bio, Aarhus, Denmark). Paired-end 
sequences from the four samples were used to construct the de 
novo C. colocynthis transcriptome assembly with default param- 
eters. Assembly reads were also assembled against the watermelon 
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Table 2. Summary details of sequences 


produced after assembly of Citrullus colocynthis reads. 




Length/Number of reads 


N50 


1 870 bp 


Average length 


1350 bp 


Minimum length 


201 bp 


Maximum length 


1956 bp 


Total number of contigs 


20581 


Full length cDNAs 


5038 


doi:10.1371/journal.pone.0104657.t002 



(C. lanatus) genome sequences, which were downloaded from the 
Cucurbit Database (http://www.icugi.org/cgi-bin/ICuGI/index. 
cgi). Reads were filtered and assembled using the CLC genomics 
workbench. The parameters used were as follows: 2 points of 
mismatch cost, 2 points of insertion cost, 2 points of deletion cost, 
0.5 as length fraction, 0.95 as similarity fraction. After the de novo 
assembly and watermelon mapping assembly, we used Trinity to 
assemble all the contigs with the default parameters. 

Identification of Full Length cDNAs 

Two methods were used to identify full length cDNAs. First, 
BLASTX searching (E value: le~'") was used to detect the 
matched cDNAs in SwissProt database; second, ESTScan 2.0 was 
used to identify the translated sequences. Sequences with either 
start codon (ATG) and stop codon (TAG/TGA/TAA), or 
sequences with start codon (ATG) and homologue to a known 
protein with S80% similarity, were chosen as full length cDNAs. 
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Figure 1. Principal component analysis of the leaf transcrip- 
tome of Citrullus colocynthis at 0 (Dl), 24 (D2), 36 (D3) and 48 
(D4) hours of withholding water. PCA analysis was conducted using 
the CLC workbench. The X-axis indicates principal component 1, the Y- 
axis principal component 2. Data from Dl and D2 are more closely 
related than data from D3 and D4. 
doi:1 0.1 371/journal.pone.01 04657.g001 



Expression Analysis Using Transcriptome Reference 

Pair-end sequencing reads of the four libraries were filtered 
using CLC Genomics workbench (0.05) before mapping to the 
references sequences from assembled cDNAs. First, read counts of 
each unigene were converted to reads per million (RPM). 
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Figure 2. Heat map of the leaf transcriptome of Citrullus 
colocynthis depicting changes in transcript patterns of 2545 
gene transcript clusters under water deficit treatment. Heatmap 
was produced using R. Each red or green block represents normalized 
expression value of a gene detected at 0 (Dl), 24 (D2), 36 (D3) and 48 
(D4) hours of withholding water. Red indicates higher gene expression 
values across treatment, while green indicates lower expression values 
across treatment. 

doi:1 0.1 371/journal.pone.01 04657.g002 
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Figure 3. Gene Ontology (GO) of significantly changed Citrullus colocynthis genes under drought stress. Based on high-score BLASTX 
matches in the NR (non-redundant) plant protein database, C. colocynthis genes were classified into three main GO categories. The y-axis indicates 
the percentage of a specific category of genes in each main category. 
doi:1 0.1 371/journal.pone.01 04657.g003 



Secondly, statistical analysis using Kal's test [22] provided in CLC 
Genomics workbench (P<0.05 and fold change ^1.5) was 
conducted. These transcripts were annotated against the reference 
sequences. 

Gene Ontology Analysis 

The functional annotation software BLAST2GO (http://www. 
blast2go.com/b2ghome) was used to conduct gene ontology (GO) 
analysis of C. colocynthis genes in this study. The databases used 
were SwissProt and NCBI. BLAST E-value was set at l.Oe ^ . The 
major GO analysis was determined by BLAST, mapping, and 
annotation. Results are presented as a bar chart showing the 
percent of genes belonging to each group. 



qRT-PCR Analysis 

For cDNA synthesis, 500 ng of the total RNA for each sample 
(the same RNA was used for RNA-seq analysis) was used in 
reverse transcription with ProtoScript First Strand cDNA synthesis 
kit from New England BioLabs (Ipswich, MA). qRT-PCR was 
performed with SYBR-Green Supermix from Bio-Rad (Hercules, 
CA) in an Eppendorf Mastercycler ep realplex (Hauppauge, NY). 
Table SI contains gene specific primer sequences. Each reaction 
contained 10 |J,1 of SYBR-Green supermix, 1 |J,1 of cDNA 
template, 1 \l\ forward primer (4 |xm), 1 ^,1 reverse primer 
(4 |j.m), and 7 |J.l ddH20. The qRT-PCR program consisted of 
one cycle at 95°C for 15 sec, followed 40 cycles of 15 sec at 95°C, 
15 sec at 55°C, and 30 sec at 72°C. The relative expression data 
was compared with actin [13] from C. colocynthis. Quantification 



Table 3. Validation of the RNA-Seq expression profiles of selected C. colocynthis genes by qRT-PCR. 



Gene ID 


Hit ID 


Annotations 


Log2FC*-RNA seq 


Log2 FC-qRT-PCR 


Compl4675 


Cla000300 


Heat shock protein 


7.60 


8.37 


Compl3927 


Cla001351 


Cold shock protein 


3.41 


4.8 


Comp2553 


Cla002576 


MA3 domain-containing protein 


2.15 


5.33 


Comp372 


Cla002814 


Bax inhibitor-1 


1.90 


3.56 


Comp8117 


Cla012143 


Drought stress related gene 


2.87 


4.31 


Compl0586 


Cla013687 


MYB 


4.09 


5.86 


Compl9751 


Cla022362 


WRKY 


6.27 


6.00 


Comp6823 


Cla022169 


MADS 


2.07 


3.67 



Transcripts identified as D4 significantly changed genes were compared to transcripts at Dl . The table shows the log2 fold change calculated from D4 expression vs Dl 
expression for RNAseq and qRT-PCR analysis. FC = fold change. Gene ID is gene named in C. colocynthis, and Hit ID Is blastx of gene with watermelon ID number. 
doi:l 0.1 371/journal.pone.Ol 04657.t003 
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Figure 4. qRT-PCR analysis of 8 selected Citrullus colocynthis genes under drought treatment. The y-axis indicates relative expression 
compared with Day 1 expression of each gene. The x-axis shows days of water-withholding time points. 
doi:1 0.1 371/journal.pone.0104657.g004 
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Figure 5. qRT-PCR analysis of 8 additional Citrullus colocynthis genes under drought treatment. The y-axis indicates relative expression 

compared with Day 1 expression of each gene. The x-axis shows days of water-withholding time points. 
doi:1 0.1 371/journal.pone.01 04657.g005 
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Table 4. List and number of several drought stress signaling pathway members detected 


in the Citrullus colocynthis transcriptome. 




Ion Signaling 


NO. 


Calcium Binding Protein 


10 


Protein Kinase Pathways For Osmotic Signaling 


MAP kinase 


16 


Osmotic Stress-Activated Phospholipid Signaling 


PLC-like Phosphodiesterase 


9 


DAG 


8 


PIP2-like Aquaporin 


5 


Detoxification Signaling 


Heavy Metal Transport/Detoxification Protein 


32 


Transcription Regulators 


Regulatory Proteins 


MYB 


131 


MYC 


43 


NAC 


42 


Leucine-rich Repeat Proteins 


150 


Functional Proteins 


Heat Shock Protein 70 


70 


Heat Shock Protein 22 1 


grpE Like Protein 


4 


RBOHD {respiratory burst oxidase) 


2 


VIP2 {VIRE2-INTERACTING PR0TEIN2) 1 


ABATransporter-like Protein 


100 


Synaptobrevin-related protein 


5 


Toc34-l (Translocon outer envelope of chloroplast) 1 


Beta-amylase 


9 


Puruvate Kinase 


11 


TIPl (TIP GROWTH DEFECTIVE 1) 1 


RD22 


2 



doi:l 0.1 371 /journal.pone.Ol 04657.t004 



of the relative transcript levels was performed using the 
comparative C-j- method. The induction ratio (IR) was calculated 
as recommended by the manufacturer and corresponds to 

2 , where AACT = (C-p, target gem-?; "Ct. a<:tm)treatmeiit" 

(Ct, target-Cx, actm)controi- All cxpcrimcnts wcrc replicated three 
times. 

Results and Discussion 

Transcriptome Assembly Results 

The transcriptome of C. colocynthis leaves following 4 days of 
drought stress was assembled and assessed following paired-end 
(2*50 bp) lUumina sequencing. The lUumina platform yielded an 
average of 24 million high-quality reads per sample (Table 1). AU 
sample reads were used to construct a de novo assembly. A 
reference assembly was constructed using the completely 
sequenced watermelon genome. A total of 20,581 contigs were 
generated (Table 2). The contigs had an average length of 
1350 bp and N50 of 1870 bp. BLASTX was used against 
SwissProt and ESTscan of translated protein sequences to detect 
full length cDNAs. A total of 5,038 full length cDNAs were 
detected in our sequencing assembly (Table 2). 



Principle component analysis (PCA) was implemented in the 
CLC workbench. The results (Figure 1) illustrate differential gene 
expression patterns in C. colocynthis seedlings following four days 
of withholding water. Gene expression patterns on Dl and D2 
were nearly similar, while results from D3 and D4 showed 
substantial differences as compared to Dl and D2, and gene 
expression patterns detected on D4 were very different from all 
other time points as a result of drought stress. 

Differentially Expressed Genes Involved in Response to 
Drought Stress in C. colocynthis 

The paired-end reads were mapped to the reference genomes 
after filtration. Read counts of each unigene were converted to 
reads per million (RPM). The read number of each cDNA was 
divided by the total number of reads per day (1, 2, 3, or 4) from the 
data set, and multiplied by 10'\ Statistical analysis was conducted 
using Kal's test (p<0.05 and fold change &1.5). Genes showing 
non-significant and significant changes in read counts are shown in 
Table S2. Each sample was compared to the day 1 sample reads 
for analysis of their significance level of gene expression. The read 
results indicated that 59 genes showed significant changes at D2, 
D3, and D4 as compared to Dl; 13 genes showed significant 
differential expression at D2 and D3 as compared to Dl; 30 genes 
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Table 5. Overview of Citrullus colocynthis genes involved in phytohormone signaling. 





Function 


Ethylene Pathway Major Members 


Number 


Ethylene-insensitive protein 


EIN 


5 


Ethylene receptors 


ETR 


3 


Constitutive triple response proteins 


CTR 


1 


EIN3-like (EIL) transcription factors 


EIL 


1 


Ethylene response factors 


ERF 


118 


Function 


Auxin Pathway Major Members 


Number 


Receptors/F-box proteins 


TIP 


1 


Ubiquitin ligase component 


SCF 


5 


Target proteins 


Aux/IAA 


30 


Auxin response factors 


ARF 


30 


Auxin transport protein 


PIN 


6 


Function 


JA Pathway Major Members 


Number 


Receptor/F-box proteins 


con 


1 


Target proteins 


JAZ 


1 


JAZ interacted proteins 


NINJA 


4 


Activator transcription factors 


R2R3-MYB transcription factor 


3 


Activator transcription factors 


MYC2,3,4 


9 


Function 


SA Pathway Major Members 


Number 


Regulatory proteins 


NPRl 


2 


SA mainly induced genes 


WRKY 


46 


SA mainly induced genes 


TGA2,3,5,6 


6 


Function 


ABA Pathway Major Members 


Number 


ABA receptors 


PYR/PYL/RCAR 


19 


PYR/PYL/RCAR interacted proteins 


PP2C 


3 


Serine/threonine-protein kinase 


SnRK2 


1 


SnRK2 target 


ABI5 


1 


Function 


GA Pathway Major Members 


Number 


GA receptors 


GID1A/B/C/-Iike 


11 


E3 ubiquitin ligase 


SLYl/SNZ 


2 


DELIA proteins 


GAl, RGA, RGLl, 2, 3 


4 



doi:l 0.1 371 /journal.pone.Ol 04657.t005 



showed significant differential expression early, at D2, as 
compared to Dl; 897 genes showed changes at D3 and D4 as 
compared to DI; 341 genes were only differentially expressed at 
D3; II91 genes were regulated late, at D4 only, under drought 
stress. In conclusion, the C. colocynthis gene expression patterns 
showed dramatic changes with 2545 genes showing significant 
changes, mostly occurring late under drought conditions (D3 and 
D4 of withholding water). 

The heat map depicted in Figure 2 corresponded to the 
principal component analysis. D3 and D4 transcripts were 
clustered together, and Dl and D2 transcripts were clustered. 
Also significant changes were seen at D3 and D4 suggesting that 
the transcriptional response of many genes was up-regulated 
during drought. Strong effects were especially observed on D3 and 
D4. 

Gene Ontology (GO) Classification 

To functionally categorize significantly changed genes in C. 
colocynthis under drought treatment, gene ontology analysis by 
BLAST2GO was performed. C. colocynthis unigenes were 
categorized in three main GO categorizes: biological process 



(2672), molecular function (1368) and cellular component (1053). 
These GO terms were further divided into several sub-categories 
(Figure 3). In the biological process category, single organism 
process genes accounted for more than 20% of the biological 
process genes. In the molecular function category, more than 40% 
of genes were associated with a catalytic activity. In the cellular 
category, more than 35% of the genes were associated with the 
cellular component. 

Validation of lllumina Expression Patterns by qRT-PCR 
Analysis 

To confirm the reliability of the lllumina sequencing read 
analysis, 8 candidate genes were selected and their expression was 
compared at D4 and D 1 using qRT-PCR. The expression patterns 
resulting from qRT-PCR showed general agreement with those 
from the lllumina sequencing analysis (Table 3). Discrepancies 
with respect to ratio of fold changes between sequencing and qRT- 
PCR analysis can be attributed to the essentially different 
algorithm and sensitivity of the two techniques [23]. In the 
deep-sequencing method the absolute expression rather than 
relative expression as in qRT-PCR analysis is used. Transcrip- 
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Table 6. Photosynthesis-related and chlorophyll-related gene expression profiles. 



Gene ID 


Hit ID 


Annotations 


D1RPM 


D2RPM 


D3RPM 


D4RPM 


Comp808 


Cla002545 


Photosystem 1 psaA 


1180 


1066.9 


513.6 J, 


418.8 i 


Comp2604 


Cla021635 


Photosystem 1 subunit II 


309 


319 


114i 


22.8 i 


Comp9446 


Cla002576 


Photosystem 1 subunit III 


470 


428 I 


97 i 


39.1 I 


Comp4444 


Cla012670 


Photosystem 1 subunit IV 


274 


219 i 


70.4 I 


8.3 i 


Comp4514 


Cla007871 


Photosystem 1 subunit IV A 


104.9 


98.1 I 


41.2 i 


9.9 i 


Comp2126 


Cla004483 


Photosystem 1 subunit V 


292.9 


281.8 i 


95.6 i 


23.5 i 


Comp4727 


Cla007940 


Photosystem 1 subunit XI 


257.2 


229.5 I 


61.4 I 


9.6 i 


Comp7482 


Cla009814 


Photosystem 1 subunit X 


225.3 


194.1 I 


60.3 I 


12.3 i 


Comp5307 


Cla011174 


Photosystem 1 subunit X 


199 


176i 


44i 


4.6 i 


Compl999 


Cla005420 


Photosystem II polypeptide 


225 


194.1 I 


60.3 i 


12.3 i 


Comp2011 


Cla008554 


Photosystem II 5 kDa protein 


138.9 


118.7 i 


45.4 i 


2.5 i 


comp7744 


Cla013942 


Photosystem II Protein 


495.8 


515.8 


331.4 i 


116.7 i 


Comp8614 


Cla014815 


Photosystem II reaction center W protein 


310.2 


265.5 


117.3 i 


28 I 


Comp2016 


Cla022723 


Photosystem II core complex proteins psbY 


78 


76.3 I 


24.7 i 


7.7 i 


Compl2757 


Cla011748 


Chlorophyll a-b binding protein 13 


290.5 


246.2 i 


3.1 i 


0.3 i 


Comp6677 


Cla013483 


Chlorophyll a-b binding protein 3C 


130.6 


89.5 I 


7.9 i 


3.8 i 


Comp1476 


Cla013826 


Chlorophyll a-b binding protein 


671.4 


583.9 I 


42.9 I 


8.8 i 


Comp7589 


Cla015680 


Chlorophyll a-b binding protein 37 


737.2 


489.3 i 


16.6 i 


10.5 i 


Comp2286 


Cla017325 


Chlorophyll a-b binding protein 3C 


577.4 


230.8 I 


4.1 i 


U 


Comp544 


Cla017983 


Chlorophyll a-b binding protein 5 


235.2 


219.1 I 


72.3 I 


48.2 i 


Comp3797 


Cla018117 


Chlorophyll a-b binding protein 6 


632.9 


502.3 I 


1 54.8 J, 


38.3 J, 


Compi 5099 


Cla019595 


Chlorophyll a-b binding protein 21 


228.8 


126.3 I 


0.1 i 


1 


Compl922 


Cla001764 


Chlorophyll a-b binding protein 8 


609.1 


470 I 


241.6 J, 


58.6 i 


Compl584 


Cla012368 


Chlorophyll a-b binding protein 8 


1082.2 


884.8 i 


111.1 i 


281.1 i 


Comp940 


Cla009752 


Chlorophyll a-b binding protein 21 


652.6 


514.8 I 


43.3 i 


3.6 J, 


Compl3569 


Cla009753 


Chlorophyll a-b binding protein 21 


827.5 


676.3 I 


68.1 1 


5.8 i 


Comp4958 


Cla022963 


Chlorophyll a-b binding protein 7 


538.2 


404.7 i 


110.1 i 


18.8 i 


Comp4043 


Cla001790 


Oxygen-evolving enhancer protein 1 of 
photosystem II 


1088.6 


795.9 I 


334.6 i 


82.8 i 


Comp3077 


Cla005429 


Oxygen-evolving enhancer protein 2, 
chloroplastic 


604.2 


531.1 I 


395.2 J, 


122.6 i 


Comp5901 


Cla019423 


Oxygen-evolving enhancer protein 3 


545.6 


414.8 i 


86.7 i 


17 i 



D1RPM-D4RMP are gene expression reads per millions for each day. I indicates the down-regulation of each gene at each time point. 
doi:l 0.1 371 /journal.pone.01 04657.t006 



tional qRT-PCR analyses of 16 genes during the drought 
treatments are shown in Figures 4 and 5. 

Gene Comp5873 is homeobox-leucine zipper protein, with 
significant upregulation during all days of drought in C. 
colocynthis. It is known that the expression of several homeobox- 
leucine zipper proteins is correlated to stress. Athb-12, a 
homeobox-leucine zipper domain protein from Arahidopsis , is 
functionally involved in salt tolerance in yeast [24]. Hahb-4, a 
homeobox-leucine zipper gene is potentially involved in water 
stress in sunflower [25]. Comp862 belongs to the glutathione S- 
transferase (GST) family, which contains heterogeneous, multi- 
functional dimeric proteins. This gene is highly up-regulated (60x) 
during drought in C. colocynthis. It is thought that GSTs are 
involved in cellular detoxification [26]. Compl 108 is a member of 
the NAG (NAM/ATAFl, 2/CUG2) gene family and NACs are 
known to be involved in numerous biological processes, including 
drought stress [27]. 



Compl0156 is a member of the GIDl (GIBBERELLIN 
INSENSITIVE DWARF 1) family in C. colocynthis with high 
expression under drought conditions. GIDl is a soluble GA 
receptor in rice [28]. GA-GIDl complex interacts with DELLA 
proteins, which are negative regulators of GA signaling pathway. 
RD22, a known dehydration responsive gene in Arahidopsis, 
which is mediated by abscisic acid (ABA), may have physiological 
and molecular significance for processes underlying memory 
functions of plants in response to ABA and light pulses [29,30]. 
One RD22-like protein from soybean can alleviate salinity and 
osmotic stress [31]. RD22-like genes (Comp 6528) with significant 
up-regulation (>160x) under drought in C. colocynthis were 
confirmed by qRT-PCR, especially after days of withholding 
water. 

NDRl (NON-RACE-SPECIFIC DISEASE RESISTANCE 1) 
plays an important role in coordinating broader ceUular processes 
in response to stress and bacterial pathogen infection [32]. The 
chemical P-aminobutyric acid, which is known to induce resistance 
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Table 7. Expression profiles of some citrulline metabolic genes during drought stress. 



Gene ID 


Hit ID 


Annotations 


D1RPM 


D2RPM 


D3RPM 


D4RPM 


Compl556 


Cla006970 


Carbamoyl-phosphate synthetase 


64.6 


70.9 t 


159 1 


175.8 t 


Comp1556 


Cla005591 


Carbamoyl-phosphiate synthetase 


196.1 


162.5 


846.5 t 


1542 t 


Comp3843 


ClaO 14787 


Carbamoyl-phosphate synthetase 


15 


20.5 t 


49.8 t 


63.3 t 


Comp 11009 


ClaO 16474 


Proline dehydrogenase 


2.2 


3.9 t 


11.6 t 


31.4 t 


Com pO 1186 


Cla023055 


Arinosuccinase 


5.6 


6.2 


3.7 


5.3 


Comp5521 


Cla020781 


Orinithine carbamoyltransfrase 


25.2 


19.1 I 


52.5 t 


32.5 


Comp6244 


Cla015337 


Acetylornithine aminotransferase 


39 


45.3 t 


127 t 


147.6 t 


Comp6587 


Cla008748 


Glutamine amidotransferase 


25.6 


19.7 


12.4 


13.8 


Compl5261 


Cla017928 


Glutamate 5-kinase 


4.3 


5.7 


1.7 


1.6 


Comp3961 


Cla019569 


Orinithine-oxo-acid transaminase 


19.3 


30.5 


28.1 


17.1 


Comp7707 


Cla023055 


Argininosuccinate lyase 


5.6 


6.2 


3.7 


5.3 


Comp3468 


Cla003592 


Argininosuccinate lyase 


22.9 


26.8 


35.4 


25.6 


Comp776_seq2 


Cla00261 1 


Arginosuccinate synthase 


12.5 


25.4 


6 


7.7 


Comp776_seql 


Cla002609 


Arginosuccinate synthase 


18.2 


18.9 


2.4 i 


2i 


Compl556 


Cla022915 


Carbamoyl-phosphate synthetase 


115.4 


103.2 


53.4 i 


42 I 


D1RPM-D4RMP 


ndicate gene expression 


reads per millions for each day. ^ and i 


indicate the significant 


up-regulation and down-regulation of each gene at each 



time point. 
doi:1 0.1 371 /journal.pone.Ol 04657.t007 



in plants, primed the expression of many genes, and NDRl/ 
NHLIO was one of them [33]. Comp73I7 gene, which is a 
member of NDRI/NHIO showed significant changes at D2 only 
during the early stage of water deficit stress. 

GRAS (for GA Insensitive, REPRESSOR of gal-3 (RGA), 
SCARECROW (SCR)) transcription factors, have a major 
function in plant development and environmental adaption. These 
TFs are particularly implicated in the modulation of plant 
tolerance to stressors as cold, drought, salinity by crosstalks via 
GA to ABA-dependent and ABA-independent pathways [34] . For 
example, SCL7 confers salt and drought tolerance in Aj'ahidopus 
[35]. SCL14 is involved in the detoxification of xenobiotics and 
possibly endogenous harmful metabolites [36]. Comp20554, 
which belongs to the GRAS transcription factor family, showed 
significant up-regulation especially on D2. 

The expression profile of Comp3048, which codes for a 
methyltransferase, maintained high levels during drought stress, 
especially on D2. It was found that myo-inositol-O-methyltrans- 
ferase (Imtl) responded to low temperature stress in transgenic 
Arabidopsis [37]. Trithorax-like H3K4 methyltransferase from 
barley is drought inducible [38] . The methylation of myo-inositol 
catalyzed by myo-inositol methyltransferase (IMT) occurs when 
plants are under abiotic stress. Over-expression of IMT resulted in 
improved tolerance to dehydration and salt stress treatment in 
Arabidopsis [39]. 

Heat-shock proteins (HSPs) are environmentally induced 
proteins that enable plants to cope with heat and other 
environmental stresses. For example, Trichoderma harzianum 
Hsp70 transgenic Arabidopsis is abiotic stress tolerant [40]. 
Similarly HSP22 was found to be highly upregulated in C. 
colocynthis roots during drought conditions [13]. Overexpression 
of GmHsp90S can decrease damage of abiotic stresses in 
Arabidopsis [41]. Compl4675 belongs to the HSP family, and 
showed up-regulation at later stages of D3 and D4 (Figure 5). 

Plant cold shock proteins (CSP) are very conserved among 
various plant genera [42]. The first CSP identified was WCSPl 
from winter wheat, which did accumulate in response to low- 



temperature stress [43]. Similarly in C. colocynthis CSP 
(Comp 13927) was up-regulated during later stages of water deficit. 

Comp2553, one MA3 domain-containing protein gene, showed 
marked changes on D3 and D4 under drought conditions 
(Figure 5). Loss-of-function of ECIPl (one MAS domain-contain- 
ing protein) resulted in enhanced ethylene response but altered salt 
response [44] in Arabidopsis. 

Overexpression of plant BL-1 in Arabidopsis resulted in the 
attenuation of cell death induced by biotic stresses (pathogens) and 
abiotic stresses such as heat, cold, drought, salt and chemical- 
induced oxidative stresses [45] . BL- 1 might function to control the 
level of the "pro-survival and pro-death" signals under multiple 
stress conditions in plants [46]. For example, cucumber BAX 
inhibitor- 1 is a conserved cell death suppressor induced by cold 
stress and a negative regulator of programmed cell death (PCD) 
[47]. Comp372 encodes one BL-1 gene, which was up-regulated 
by drought at later stages (D3 and D4). 

Comp8117, which is a homologue of a Populus EST 
(CU233481.1), is a drought stress related gene, up-regulated at 
D4. Comp 10586, which encodes a member of the MYB 
transcription factor family, is up-regulated during the late stage 
of drought stress. Some MYB members have been shown to 
regulate plant responses to biotic and abiotic stress conditions [48] . 
For example, AtMYB96 acts through the ABA signaling pathway 
to induce pathogen resistance by promoting salicylic acid 
biosynthesis, and thus regulating stomata movement, drought 
tolerance and disease resistance in Arabidopsis [49,50]. MYB88 
might function directly or indirectiy, as positive regulator of stress- 
responsive genes [51]. TaMYB30-B from wheat did improve 
drought stress tolerance in transgenic Arabidopsis [52]. Another 
family of transcription factors, the WRKY family, named for the 
WRKY domain of about 60 amino acids, contains a highly 
conserved WRKYGQK heptapeptide at its N-terminus and a 
zinc-finger-like motif at its C-terminus [53]. WRKY transcription 
factors are involved in multiple aspects of plant growth, 
development and stress [54]. Several TaWRKY in wheat with 
roles in the abiotic stress response acted in an ABA-dependent 
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manner [55]. Here, gene Compl9751encodes a WRKY gene, 
which showed remarkable up-regulation at D4. 

MADS-box family members function in reproductive develop- 
ment and stress [56]. For example, OsMADS25 and OsMADS27 
transcripts accumulate in response to osmotic stress [57]. 
Comp6823, which encodes for a MADS-box gene in C. 
colocynthis, is up-regulated markedly at the last stage of drought 
(D4; Figure 5). 

Analysis of the Drought Stress Signaling Transcriptome in 
C. colocynthis 

Drought stress signal transduction consists of several pathways 
including ionic and osmotic homeostasis signaling, detoxification 
response and growth regulation pathways [58] . Genes detected in 
C. colocynthis leaves during water deficit are listed in Table 4. 

In the ion signaling pathway, calcium binding proteins are well 
known for their involvement in both biotic and abiotic stress 
response pathways [59]. The calcium ion (Ca^"^ as a secondary 
messenger in plants is sensed by calmodulins (CaMs)/CaM-like 
protein (CMLs), the caldineurin B-like proteins (CBLs) and Ca2+- 
dependent protein kinases (CDPKs). CaM binds to CaM-binding 
proteins (CBPs), which function in different pathways under biotic 
and biotic stresses [60]. A total of 10 Ca^"*" binding proteins were 
detected in the C. colocynlhis transcriptome. 

Protein phosphorylation is a central theme in the cell's response 
to stress. The MAP kinase cascade in transcript levels consist of a 
number of protein kinases, such as two-component histidine 
kinase, MAPKKK, MAPKK, MAPK etc. [61]. Here we detected 
16 MAP kinases in the C. colocynthis transcriptome. 

Membrane phospholipids can activate several types of phos- 
pholipases that cleave certain phospholipids to generate lipid 
messengers (eg. PA, DAG, IPS), which further regulate stress 
tolerance through modulation of stress-responsive gene expression 
[62]. Several members in this pathway, such as phospholipase C 
(PLC), diacylglycerol (DAG) and phosphotidylinositol 4,5-bispho- 
sphate (PIP2)-like aquaporin were detected (Table 4). 

Detoxification signaling can ameliorate the damage in plants 
under stresses [63], as noted in many other plant species. A total of 
32 detoxification proteins were detected in the C. colocynthis 
transcriptome. 

Molecular mechanisms regulating gene expression in response 
to drought stress have been studied by analyzing the functional 
transcription factors in ABA-dependent and ABA-independent 
pathways [6]. Several regulatory proteins in ABA-dependent or - 
independent pathway were detected, among which NAG, MYB/C 
and leucine-rich repeat proteins (LRR). These regulator^' proteins 
can further modulate many responsive transcription factors. 
Several functional proteins such as heat shock protein (HSP) 70, 
HSP22, grpE like protein, respiratory burst oxidase D (RBOHD), 
VIRE2-Interacting protein2 (VIP2), ABA transporter-like protein, 
synaptobrevin-related protein, translocon outer envelope of 
chloroplast (Toc34-l), beta-amylase, TIPl (TIP GROWTH 
DEFECTTVE 1), and RD22 were detected. 

Analysis of Phytohormone Signaling Mediators in 

C. colocynthis 

Phytohormones play important roles in regulating plant 
responses under biotic and aljiotic strc-ss. Elaborate phytohormone 
signaling networks mediate the adaptability of plants to difierent 
environmental conditions [64]. Many phytohormones such as 
ABA, salicylic acid (SA), jasmonic acid (JA), auxin, ethylene and 
gibberellic acid (GA) are being studied for their role in abiotic 
stress responses [33,65,66]. It is known that downstream signaling 



proteins for auxin, GA, JA, and ABA are subjected to ubiquitin- 
dependent degradation [67]. Putative phytohormone signaling 
genes detected in C. colocynthis during the drought response are 
listed in Table 5. 

Ethylene signaling pathway components were ordered into a 
hypothetical linear pathway based on both genetic (epistasis) 
analysis and biochemical interactions [68]. Almost all of the 
ethylene signcding homologous members (118 ethylene response 
factors or ERFs) were detected in the C. colocynthis transcriptome. 
APETALA2/ ethylene responsive factor (AP2/ERF) transcription 
factors are well-known for mediating stress responses and 
development in plants [69] . 

The auxin response factor (ARE) family contains transcription 
factors that bind to auxin-responsive elements (AREs) in the 
promoters of primary auxin-responsive genes. Aux/IAAs are early 
auxin-response proteins that bind ARFs, therefore inhibiting 
ARE-mediated gene transcription. Aux/IAAs are involved in 
ubiquitin-mediated degradation, which is catalyzed by SCF E3 
ubiquitin ligase. TIRl can stimulate Aux/IAA proteolysis by 
binding auxin to this protein [70]. All of the major auxin signaling 
related transcription factors found in the C. colocynthis transcrip- 
tome are shown in Table 5. Similar to the auxin pathway, a novel 
family of transcriptional regulators, the jasmonate ZIM-domain 
(JAZ) proteins play a target part as Aux/IAA. In addition, COIl 
plays a similar role as TIRl, while MYC and R2R3-MYB 
transcription factors work as ARFs [71]. Several JA signaling 
pathway related genes exist in the C. colocynthis transcriptome. 

In the SA pathway, NPRl (non-specific disease resistance 1) is a 
key regulator in SA-dependent defense signaling [72]. Similarly 
WRKY and TGA play major roles as transcriptional regulators in 
the SA patiiway. We detected 2 NPRl related proteins, 46 
WRKYs and 6 TGAs in C. colocynthis, which might function in 
the C. colocynthis SA pathway. 

Drought triggers the production of the phytohormone ABA, 
which in turn causes stomatal tiosure and induces expression of 
stress-related genes. The soluble PYR/PYL/RCAR receptors 
function at the apex of a negative regulatory pathway to directiy 
regulate PP2G phosphatases, which in turn directiy regulate 
SnRK2 kinases [73]. Several of the core transcription factors in 
the ABA pathway are listed in Table 5. 

We also detected 1 1 GA receptor GIBBERELLIN INSENSI- 
TIVE DWARF 1 (GIDl), 4 DELIA growth inhibitors (DELLAs) 
and 2 F-box proteins (SLYl) and SNEEZY (SNZ), which play 
important roles in GA signaling pathways [74,75]. 

Cellular Metabolism under Drought Stress in 
C. colocynthis 

Global gene expression analyses have shown substantial down- 
regulation of many photosynthetic genes under drought not only in 
Arahidopsis [76], but also several other species such as indica rice 
[77,78]. Similarly, many photosystem I and II, chlorophyll a, b binding 
protein, and oxygen evolving enhancer protein genes [79], showed 
down-regulation during water deficit stress in C. colocynthis (Table 6). 

Citrulline, a non-protein amino acid intermediate in the aipnine 
biosynthetic pathway, has been found to accumulate in leaves of 
drought tolerant watermelon under water deficit conditions [80,8 1,82]. 
Thus several factors related to die response of this model drought 
tolerant species to stress have been identified (Table 7, Figure SI). 
Citrulline metabolic genes (carbamoyl-phosphate synthetase, acetyl 
glutamate synthase, acetylomithine aminotransferase, aminoglutamate 
decarboxylase, acetylomithine deacetylase, and glutamate dehydroge- 
nase) were found to be significantly up-regiilated during drought. 

One of the major research goals is to understand the molecular 
mechanisms underlying drought tolerance in plants. It is clear that 
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drought triggers a wide variety of responses in C. colocynthis. 
Down regulation of many photosynthetic genes was observed 
especially at the later stages of drought. Up regulation of many 
transcription factors, stress signaling factors, detoxification genes, 
and genes involved in phytohormone signaling occurred through- 
out the water deficit experiment. Systematic approaclies using 
genomic analyses should lead to the discovery of additional stress 
factors and provide us with a better understanding of the stress 
tolerance mechanism of this drought tolerant plant species. 

Supporting Information 

Figure SI Citnilline metabolism pathway in C. colo- 
cynthis. 

(TIF) 
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